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ABSTRACT 

We present calculations of expected continuum emissions from Sedov- Taylor 
phase Type la supernova remnants (SNRs), using the energy spectra of cos- 
mic ray (CR) electrons and protons from nonlinear diffusive shock acceleration 
(DSA) simulations. A new, general-purpose radiative process code, Cosmicp, 
was employed to calculate the radiation expected from CR electrons and pro- 
tons and their secondary products. These radio. X-ray and gamma-ray emis- 
sions are generally consistent with current observations of Type la SNRs. The 
emissions from electrons in these models dominate the radio through X-ray 
bands. Decays of tt^s from p — p collisions mostly dominate the gamma-ray 
range, although for a hot, low density ISM case {njsM = 0.003 cm^'^), the 
pion decay contribution is reduced sufficiently to reveal the inverse Compton 
contribution to TeV gamma-rays. In addition, we present simple scalings for 
the contributing emission processes to allow a crude exploration of model pa- 
rameter space, enabling these results to be used more broadly. We also discuss 
the radial surface brightness profiles expected for these model SNRs in the 
X-ray and gamma-ray bands. 
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1 INTRODUCTION 



Diffusive Shoclc Acceleration (DSA) at astr ophysical shocks has bec ome the standard model 

and references therein). 



Malkov fc Drury 



2001 



for production of cosmic rays (CRs) {e.g., 
The CRs observed at the knee (~ lO^^eV) and below are most c ommonly explained through 



DSA in the blast waves of Galactic supernova remnants (SNRs) (jBlandford fc Ostriker 



Lagage fc Cesarsky 



19831: 



Blandford &: Eichler 



1987 



Drury et al. 



1978 



20011). Several indirect ar- 



guments, iiicluding energetics a nd composition, support this explanation {e.g.. 



Gaisser 



2005 



Hillas 



2006 



Ellison et al. 



20071). 



Direct evidence that at least some SNRs are capable of accelerating electrons and possi- 



bly hadrons to at least tens of Te V comes 
some historical SNRs {e.g., iBamba et al. 



rom observations o 



2006 



Funk et al. 



Parizot et al. 



2007 



nonthermal X-ray filaments in 



20061) and from TeV 7-rays ob- 



Lemoine- Goumard 



served in several shell SNRs {e.g., 
mal X-rays are likely to be synchrotron emission from TeV electrons {e.g., 



2007h. The nonther- 



Uchiyama et al. 



20071 ). The TeV 7-ray origins, although they obviously require TeV or higher energy charged 



particles, have less certain origins and are sensitive to several model parameters. Inverse 
Compton scattering, secondary 7r° decays and bremsstrahlung may all contribute, but at 
levels that depend on details of the accelerated CR spatial and energy distributions as well 
as ambient photon, magnetic field and plasma properties. However, recent observations by 
Fermi of several Type II SNRs indicate that secondary vr^ decay from proton-proton interac- 
tions is the most likely explanation for th e 7-ray emission t here, which would observationally 
confirm the SNR origins of CR protons (lAbdo et al.l 120101 ). 

In order to test the SNR origins of Galactic CRs, generally, and the accuracy and 
physical assumptions of DSA model calculations, specifically, it is important to determine 
from simulations not only the CR distributions, but also the relate d emissions, so they 



can be compared with current and anticipated observations {e.g., 



198 



3 



Blandford fc Eichler 



2009 



Villas 



20051). Indeed several model SNR comparisons of this type have been published 



Ellison &: Cassam-Chenai 



Zirakashvili fc Aharonian 



2005 



Berezhko. Ksenofontov fc Volk 



2009 



Morlino. Amato fc Blasi 



2OIOI ). Those demonstrate that observed emissions from 
some individual SNRs can be explained by the simulations with reasonable model choices, 
although they do not yet provide unambiguous model tests nor clear confirmation of the 
specific origins of TeV emissions. 

At this early stage of our understanding it is especially useful to explore generally the 
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nonthermal emissions produced by CRs in common classes of SNR models and how they 
depend on model parameters. In this spirit we present here an examination of the nonthermal 
radio to 7-ray emissions for Sedov- Taylor phase Type la SNRs, using a set of time dependent 
simulations that followed the evolution of spherical blasts, including nonlinear DSA at quasi- 
parallel shocks with thermal leakage injection of CR protons and electrons out of the bulk 
plasma. The pressure of CR protons is included in the gasdynamic equations, while CR 
electrons are treated as a test-particle component. The principal differences expected in 
the electron and proton CR populations come from differences in injection rates from the 
thermal population and from significant synchrotron/Compton radiative losses experienced 
by the highest energy electrons, but not the protons. We consider four different SNR models, 
including a range of uniform external environments and blast energies. 

The structure of our paper is as follows. Section 2 provides a brief outline of the simula- 
tions that include both CR protons and electrons. A description of the code used to calculate 
nonthermal emissions as well as results from our emission calculations are presented in §3, 
while ^4 summarizes our conclusions. 



2 NONLINEAR DSA SIMULATIONS OF SN lA SEDOV-TAYLOR 
REMNANTS 

2.1 The Spherical CRASH code 

In the simulations discussed here the evolution of the CR modified SNR shock is followed 
using a one-dimensional, spherically symmetric version of the Cosmic RAy SHock (CRASH) 
gasdynamic-CR code. The basic gasdyna mic equations modifie d to implement nonlinear 



DSA and full code details can be found in 



Kang fc Joned (120061 ). CRASH applies adaptive 



mesh refinement (AMR) techniques and subgrid shock tracking to obtain high spatial reso- 
lution close to the blast shock, where it is crucial for converged solution of the coupled gas 
dynamics and CR transport equations. The spherical CRASH code also incorporates comov- 
ing coordinates expanding with the SNR blast. The momentum- dependent CR distribution 
is evolved through a time-dependent solution of the diffusion convection equation (DCE). 
CRASH incorporates dynamical backreaction from CRs onto the gas dynamics, including 
influence of the CR pressure, Pc-, and energy exchanges coming from injection of low en- 
ergy CRs at the gas subshock and also dissipation of Alfven waves stimulated by streaming 
CRs. We note that Berezhko and collaborators have extensively studied the similar problem. 
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implementing a different, unconventional numerical method for the gasdynamics that nor- 
malizes the spatial variable by diffusion length at each momentum value and solves the CR 
kinetic equat ion iteratively to match the downstream and upstream solutions at the sub- 
shock {e.g., 



Berezhko et al.l 



1994 



Berezhko fc Volk 



1997; 



Berezhko. Ksenofontov fc Volk 



20021 ). We find that the results of our simulations using the more conventional CRASH code 



are quite consistent with these previous studies. 

Since DSA operates on relativistic electrons of a given p in exactly the same way as it 
does on protons (i.e., it is exclusively rigidity dependent), the pitch-angle-averaged phase 
space distribution functions for CR proton and electron components, fp{p, r, t) and /e(p, t, t) 
both obey the DCE, 



2 / ^^9 
r K{r,y) — 



dg , .dg 1 d r n. .1 fdg \ Id 
dt dr Sr"^ dr ^ \dy J r^ dr 

Here g represents either gp = p^fp or = P^/p, V = ln(p) and K(r, y) is the spatial diffusion 

coefficient in the radial direction (jSkillinglll975l ). Henceforth, particle momenta, p, of both 

protons and electrons are expressed in units of m^c, with the proton mass. 

The cooling term ish{p) = —dp/dt = for protons and 

rrip (TtC ^ 



Hp) 



nip Gvr nipC 



(2) 



-2.37 X 10~^^s"^ 



cS 



p 



100 /iC, 

for electrons combining synchrotron and IC cooling, where cr^ is the Thomson cross section, 
and -Bpg = BirY + B^ is the effective magnetic field strength including the energy density 
of the ambient radiation field. We discuss the magnetic and radiation energy densities in 
§3.2. The code does not include direct backreaction from large scale magnetic fields, nor the 
pondermotive force of the wave turbulence. 

The velocity m„, represents the effective radial motion of scattering centers with respect 
to the bulk flow velocity, u. Assuming that waves upstream of the gas subshock are pre- 
dominantly resonantly generated through CR streaming and that the shock has a parallel 
magnetic field geometry, upstream of the gas subshock u^, is set to the Alfven speed, va> ^■ 
The current version of CRASH does not follow the self- consistent evolution of the magnetic 
field strength through wave-particle interactions. So here we simply set va = Bq/ A/47rp(r), 
where Bq is the upstream mean magnetic field strength. Downstream, the simulations as- 
sume = 0, since the Alfvenic turbulence in that region is probably relatively balanced. 
The transition in scattering center motion, Au^, > 0, reduces the effective velocity difference 
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experienced by CRs across the shock compared to the bulk flow, Au < 0, thus reducing the 

DSA efficiency. Gas heating due to Alfven wave dissipation is represented by the term 
dP 

Wir,t) = -UHVA^, (3) 



dr ' 

where uh is a commonly used dimensionless parameter that controls the degree of dissipa- 
tion. This dissipation term derives from a simple model in which Alfven waves are reso nantly 



amp 



2006 



ified by streaming CRs in balance with local w ave dissipation processes (e.g. 



1993 ) . As previously shown in SNR simu 



Caprioh et al. 



.e.g. 



201 



McKenzie fc Volk 



Ptuskin et al. 



1982 



ation s {e.g., 



Berezhko fc Volk 



1997 



Jones 



Kang fc Jones 



20101). and in more general contexts in earlier work 



Markiewicz. Drury. fc Volk 



1992 



Jones 



19931 ). Alfvenic drift 



and precursor heating by wave dissipation both reduce DSA efficiency and associated shock 
modification. The significance of Alfven wave dissipation and drift can generally be charac- 
terized in terms of the Alfvenic Mach number, i.e., va = uq/Ma, where uq is the upstream 
plasma flow speed relative to the shock. In these simulations Ma ^ 1, so the influence is 
modest. 

Bohm-like spatial CR diffusion in the radial direction was used for both protons and 
electrons in the simulations; namely. 



Po 
p{r) 



(4) 



where = mpC /(3ei?o) = 3.13 x 10 cm s^ , is the upstream magnetic field 
strength in micro-Gauss, and po is the upstream density. The density dependence in this 
diffusion model represents enhancement of resonant Alfven wave amplitudes through com- 
pression. 



Low energy CR proto ns and electrons were in 



through thermal leakage ( IGieseler. Jones, fc Kang 



ected at the shock in the same manner 



2000 



Kang. Jones, fc Gieselerll2002l ). In 



this model, thermal ions in a Maxwellian distribution immediately downstream of the gas 
shock have a finite probability to escape upstream into the low energy CR population, pro- 
vided they have sufficient upstream-directed momentum to overcome transverse MHD wave s 
generated through the cyclotron instability in the shock transition ( iMalkov fc Vokl 119981 ). 
This behavior is modeled numerically through a "transparency function" . The function con- 
tains one adjustable parameter, = Bq/B±, which compares the upstream magnetic field, 
Bq, with the amplitude of postshock, amplified MHD waves that interact with low energy 
particles, Bj_. The value = Bq/Bj_ = 0.25 was used for the models considered here. The 
resulting fraction of thermal protons injected into the CR pool, ^ ~ 10~^. It was shown pre- 
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viously that the DSA efficiency saturates at such a high injection rate (if ^ > 10~^ for strong 
shocks), so except for details near start up, the SNR simulation results are insensitive to e^. 
The temperature used for the thermal population was computed from the gas pressure, Pg, 
and density, p, as T = finipPg/ (ksp) with p = 0.61, so assuming equal proton and electron 
temperatures. 

Because postshock thermal electrons have smaller gyroradii, compared to thermal pro- 
tons, the injection rate of electrons should be much smaller. The ratio of CR electron to 
proton number injected at the shock, K^/p ~ 0.01, is commonly assumed, since about 1% of 



the o bserved Galactic CR flux near a GeV is due to electr ons (jSchlickeiser 



20081 ). On the other hand. 



2002 



Reynolds 



Morlino. Amato fc Blasil ( l2009l ) point out that the CR electron 



flux observed here at Earth is actually a convolution of CR electrons from SNR and other 
electron accelerators of all ages. It is also entirely possible that CR electron injection, which 
is not well understood, and acceleration efficiency change dramatically for SNRs of different 
ages. On the other hand, it is commonly found in m odels of young SNRs that the best flt 



to the observational data requires K, 



Morlino. Amato fc Blasi 



2009 



10" 



Berezh 



Zirakashvili fc Aharonian 



TO. Ksenofontov &: Volk 



2009; 



20101), so we set K^/p = 10 ^ in our 



calculations. However, since the electron population is passive and the SNRs are optically 
thin in emission bands of relevance here, our resulting emissions can be generally linearly 
scaled to another preferred value of Ke/p. 

In flows where CR backreaction is important, the formation of a CR precursor compresses 
and heats the inflowing plasma. These developments can lead to substantial changes in the 
strength of the dissipative, gas subshock and to the postshock conditions relative to those 
in a pure gas dynamic shock of the same Mach number. Given these flow modifications in 
front of the subshock, it is helpful in the following discussion to identify specifically the 
unmodified, upstream conditions by the subscript '0', the conditions immediately upstream 
of the gas subshock by the subscript '1', and the conditions immediately downstream of the 
full shock structure by the subscript '2'. 



2.2 Remnant model parameters 

Table 1 lists the dynamical parameters for the SNR models: the uniform ambient density, 
riisMi the ISM temperature, Tjsm, the upstream magnetic field strength, Bq, and the SN 
explosion energy Eo- The parameters for SI- S3 (S4) represent warm (hot) phase ISM en- 
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Table 1. SNR model parameters 



Model 


niSM 


TiSM 


Eo 


Bo 


ro 


to 


Mo 




(cm^'^) 


(K) 


(lO^^ergs) 




(pc) 


(years) 


(10* km/s) 


SI 


0.3 


3 X 10* 


1 


30 


3.19 


255. 


1.22 


S2 


0.3 


3 X 10* 


4 


30 


3.19 


127. 


2.45 


S3 


0.3 


3 X 10* 


1 


5 


3.19 


255. 


1.22 


S4 


0.003 


10^ 


1 


5 


14.8 


1182. 


1.22 



Note: The model ISM Alfven speed, va = 101 km/s for S1,S2; da = 16.8 km/s for S3; va = 168 km/s 
for S4. 

vironments. Models SI, S3, and S4 assume an explosion energy, Eo = lO^^erg, while model 
S2 begins with a blast containing four times this energy, so = 4 x lO^^erg. All models 
assume there are no pre-existing CRs in the ambient ISM. 

The SN ejecta mass is set to Mgj = 1.4Mq. The simulations are intended to follow 
evolution during the adiabatic, Sedov- Taylor (ST) evolution stage. Thus, convenient nor- 
malization variables are po = 2.34 x 10~'^^njsM g cm~^ for density, along with ST similarity 
scales ro = [3Mej/{A7rpo)f^ for length, and to = {porl/Eof^ = [{3M,j/{A7r)f/^po'^''E;'^^ 
for time. The velocity and pressure scales are defined as Uo = ro/to = ^/AttE^JJsM^ and 
Po = PoUo = Eo/r^, respectively. In order to avoid confusion later, we note here a distinction 
between ST normalization subscripts, 'a ' and initial or upstream condition subscripts '0 

It is worth noting that the mass swept up by the blast at time t and radius = 
is'^oit/toY^^ is given in the ST solution as 

Ms^{t) = M,,es{rslrof = M,,esit/tof'\ (5) 

where ~ 1-15 is the ST similarity constant for a blast in a uniform medium with the 
gas adiabatic index 7^ = 5/3. The ST shock speed can be expressed as 

- - {£) - (£) - «x iw^) iB 

The second to last relation reveals that Ust Eo/Mgw, emphasizing that the dynam- 
ical state of a ST blast is determined by the blast energy and the mass contained within the 
blast. Thus, except for numerical factors of order unity, many of the results we present below 
can be approximately represented independent of the ambient radial density distribution in 
terms of M<,^ in combination with Eo and Mgj rather than t/to and Uo- 

Recent X-ray observations of young SNRs reveal in many cases mag netic fields stronger 
by at least an ord er of magnitude than the average ISM field {e.g., 



Bamba et al. 



Parizot et al 



2006; 



20061 ) . The existence of such strong fields is now commonly inte rpreted as the 



result of amplification within the shock precursor, either through resonant (ILucek fc Bell 
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2OOOI ) . or nonresonant (jBel]||2004j ) streaming instability, or possibly hydrodynamic instabilit y 



driven by the CR pressure gradient in the precursor (IBeresnyak. Jones. &: LazarianI 120091 ) . 
The simulations discussed here do not follow self-consistent amplification in the precursor 
of the magnetic field strength through these processes. Instead, to provide field values in 
the model SNRs consistent with observations, the upstream magnetic field strengths in SNR 
models SI and S2 were set to the relatively large value, Bq = 30yuG. For comparison, the S3 
and S4 models adopted, Bq = 5fiG, which is similar to the mean ISM magnetic field. 

In the simulation of the electron synchrotron energy losses and the resulting synchrotron 
emission, we assume that the field strength is modified during passage through the shock 
and into the SNR interior in a way that maintains consistency with the diffusion coefficient 
model (equation |4]). This k model assumes in the relativistic limit that the scattering length, 
A oc p/{pBq), along with a fixed ratio of wave field to total field strength. Consequently, it 
assumes B oc p, so that gas compression through the modified shock, x = P2/P0 ~ 4 — 13 
(see Fig. 2), leads to postshock field values in these models, B2 = x^o ~ 120 — 400 /iG in 
SI and S2 and ~ 20 - 50 pG in S3 and S4. 

In all the simulations CR backreaction quickly increases compression from the initial 



X = 4tox~7 — 13. Subseque ntly, y decreases s. 



Kang et al. 



owly in each case roughly as x oc t 



(120091 ) found for strong, CR modified, plane 



with k 1/7 - 1/5 (see Fig. 2) 
shocks that approximately x M.^^^, where Ai is the shock sonic Mach number. For these 
SNR simulations, Ai oc t^^/^, so their result would predict k ~ 2/15, which is reasonably 



consistent. For ST SNRs one has more generally that Ai (x M. 
discussed in §3 we will assume that x ^ Msw^^ ■ 



-1/2 



so for scaling relations 



The SNR simulations start with ST similarity blast waves at a time t/to = 0.5, because 
the detailed evolution before that time does not strongly affect the later development. In 
fact, hydrodynamic simulations of SN blast waves (without the CR terms) show that the 
evolution of the outer shock speed can be approximated by Ust for t/to > 0.2, although the 
true ST solutio n is establish ed only after the inner reverse shock is reflected at the center at 
t/to ~ 7, {e.g., 



Kang 



20061 ). The simulations are carried out to t/to = 15, which is deemed 



sufficient to establish basic ST phase properties of the CR population and their dynamical 
impact until either the blast strength weakens substantially or the blast becomes radiative. 
Earlier studies showed that the highest momentum produced by DSA during the expansion 
of an SNR, Pmax, is achieved near the end of the free expansion stage, so when t/to ~ 1, and 
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Figure 1. The time evolution of gas density and temperature for model SI along with CR pressure and the factor p^T^/'^, 
which provides a scaling for thermal bremsstrahlung. Times shown are: t/to = 1 (red dotted), 3 (green short dashed), 6 (blue 
long dashed), 10 (black solid), and 15 (magenta dot-dashed). (A colour version of this figure is available in the online journal.) 

the transfer of explosion energy to the CR component oc curs mostly during the early ST 



stage, so also when tjto ~ 1 {.e.g., 



Berezhko &: Volk 



19971). 



2.3 Evolution of the CR modified SNRs 

We now summarize the distinctive time evolutionary features of the SNR models used in 
this study, focusing on nonlinear DSA influences that impact radiative emissions. Fig. 1 
shows radial profiles from model SI of several useful quantities near the shock at times 
tjto = 1, 3, 6, 10, 15. The SI model is representative in this regard of all the simulations. Gas 
density and temperature are shown in the top two panels. They are combined in the lower 
right to show p^T^/^ as a proxy for bolometric thermal bremsstrahlung emissivity. The CR 
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pressure profile is included in the lower left panel. To simplify the plots, initial ST conditions 
from t/to = 0.5 are not included. At the first time shown in Fig. 1, t/to = 1, the shock is 
modified substantially through compression and heating in the CR precursor. In fact, the 
CR acceleration and nonlinear modification peak early in the evolution and then decrease 
in time as the shock slows down (see also Fig. 2). 

CR modification of the shock enhanced the postshock density but reduced the gas pres- 
sure, leading to a postshock gas temperature lower than that expected for the ST shock 
near t/to ~ 1- The effect of CR pressure on the temperature distribution behind the shock 
can be seen by noting that the slope of temperature rise towards the interior steepens as 
Pc{r) begins to drop sharply inside the blast. This interior drop in Pc{r) is an artifact of the 
location of the contact discontinuity separating matter interior and exterior to the initial 
conditions shock. Its shape depends on the rate of spatial CR diffusion into the blast cavity. 

In addition to raising the total density jump within the shock, compression through the 
precursor also preheats gas before it enters the subshock structure. This adiabatic effect 
by itself would raise the upstream temperature by a factor Ti/Tq ~ (pi/po)^^^, or about 
1.6 in the SI model. In fact. Fig. 1 shows that in the SI model Ti/Tq ^ 10 — 100. The 
additional heating comes from MHD wave dissipation (equation [3]). Compression through 
the precursor also decelerates upstream flow before it enters the subshock. All these effects 
greatly weaken the subshock. So the sonic Mach number of the subshock decreases quickly 
to Aisub ~ 4.8 in the SI and S4 models and to Aisub ~ 3.7 in S3 model by t/to ~ 0.7 and 
remains approximately the same afterward, even though the total shock sonic Mach number 
decreases roughly as M oc {t/to)~^^^- Since Aisub stays constant throughout the evolution, 
the CR injection fraction (C,) via thermal leakage also stays constant after t/to ~ 0.7 (see 
Fig. 2). With the adopted injection parameter (e^ = 0.25), the CR injection fraction is high, 
^ ^ 10^'^, and so the CR acceleration should be close to maximum efficiency. 

These shock modifications also significantly influence the thermal bremsstrahlung emis- 
sion in the SNR, whose distribution is represented by the p^T^/^ profile in the lower right 
panel of Fig. 1. The thermal bremsstrahlung shell will appear thinner in the strongly mod- 
ified SNRs than in the analogous gas dynamic SNR. The enhanced postshock density also 
enhances the thermal bremsstrahlung emissivity. The thinner shell and higher emissivity in 
the shell combine for a relatively small difference in the bolometric thermal bremsstrahlung 
luminosity compared to an unmodified shock. The thermal bremsstrahlung luminosity in a 
fixed spectral band below the spectral cutoff should generally be enhanced, however, since 
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Figure 2. The time evolution of the compression factors, pi/po and x = P2/PO1 postshock pressures, Pg^2 and Pc, 2i the 
injection fraction, ^, and the volume integrated energy, E^f^ and Ec is shown for six models. The wave heating parameter is 
ljh = 0.5 for Sir, S3r and S4r models, while luh = 1- for SI, S3, and S4 models. The r models are the dashed lines. (A colour 
version of this figure is available in the online journal.) 

that has an approximate 1/T dependence with respect to the bolometric emission. Of course, 
the lowered temperatures behind the modified shock also move the upper cutoff of the ther- 
mal bremsstrahlung SED from the hard X-ray band into the soft X-ray band. 

More insights to the evolutionary behaviors of the SI, S3 and S4 models can be obtained 
from Fig. 2. As already noted, both the precursor compression, pi/po, and the total shock 
compression, x = fh/Po, as well as the CR pressure relative to PoUgj^ increase quickly from 
their initial values □, reach their maximum values around t/to ~ 0.7, and decrease slowly in 
time as the shock slows down. In case of the S2 model with higher Eo (not shown in Fig. 2), 



^ The very early evolution is affected by numerical start-up of the simulations. 
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the CR acceleration and the nonhnear modification are shghtly more significant compared 
to the SI model, because of a faster shock speed. If the Alfvenic drift and heating were not 
included, the temperature of the ambient medium (and so the sonic Mach number) would 
be the primary parameter that determines the CR injection and acceleration efficiencies. 
Compare, for instance, the S3 model in warm ISM to the S4 model in the hot ISM. Because 
of the ISM temperature differences, the sonic Mach number of the total shock in the S3 
model is 6 times higher than S4 at a given t/to- So the S3 SNR would be the more efficient 
CR accelerator. The low density and high temperature ISM of the S4 model produce low 
sonic and Alfvenic Mach numbers {M. oc T~^/^, AAa = Us/va nisM^^'^ / B), making the 
CR acceleration much less efficient in that model. 

The magnetic field strength also has a significant infiuence on the SNR evolution even 
though no direct MHD dynamical effects are modeled. Most important to these simulations, 
the field strength sets the CR spatial diffusion rate and thus, in concert with Uo, the CR 
acceleration rate (see §2.4, below). Consequently, particles are accelerated to several times 
higher momenta (pmax) in the SI model in comparison to the S3 model because of the 
enhanced magnetic field and smaller diffusion coefficient (see Fig. 3 below). On the other 
hand, the stronger field leads to a faster Alfven speed, which tends to reduce DSA efficiency 
by adding entropy to the precursor. This reduces the effective velocity jump across the 
subshock and also reduces the subshock Mach number. Those effects are responsible for 
the moderate reductions in CR acceleration efficiency in the SI model compared to the S3 
model. The effect of Alfven wave heating is further illustrated in Fig. 2, where we show 
three additional simulations. Sir, S3r, and S4r, in which the wave dissipation parameter is 
reduced from uh = 1 to uh = 0.5. The smaller wave dissipation rate reduces non-adiabatic 
heating in the precursor, allowing greater compression through the precursor, so an increase 
in CR acceleration efficiency. 

Fig. 2 also shows the volume integrated thermal energy {Eth = J 47rr^(ir(3/2)Pg) and CR 
energy in units of the SN explosion energy. The energy transfer to CRs seems to saturate 
for t/to ^ 10, and the fraction of the blast energy transferred to CRs (Ec/Eo) is about 0.7, 
0.5, 0.45, 0.35 in the S3, S2, SI, and S4 models, respectively. 

In summary, adopting simple models for Alfven wave transport, these simulations demon- 
strate that the CR acceleration would be less efficient at SNRs in hot rarefied ISM where 
both the sound speed and Alfven speed are faster, compared to those in a warm phase ISM. 

We note that the evolution of CR modified SNRs shown in Fig. 2 (especially, pi, p2, Pg,2 
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Figure 3. Volume integrated proton spectra in terms of particle kinetic energy at t/to = 1 (red dotted), 2 (green short 
dashed), 3 (blue long dashed), 6 (magenta dot-dashed), and 10 (black solid). (A colour version of this figure is available in the 
online journal.) 



and -Pc,2) differ somewhat from tlie similar SNR models presented in iKang ( 120061 ) . Several 
parameters were different there, e.g., /i = 1, = 0.16 — 0.2, Tjsm = lO^K for Sl-3 models, 
and Bq = 30/iG for the S4 model. The main differences result, however, from a coding error in 
those earlier simulations, which caused the Alfvenic drift and heating effects to not be fully 
implemented. The results in the present studv are, however, con sistent with independent 
simulations of similar SNRs ( iBerezhko. Ksenofontov fc Volkll2002l ). in which a Type la SNR 
with Eo = 3 X lO^^erg and ^ = 2 x 10~^ in warm phase ISM with ujsm = 0.3 cm~^, 
TisM = lO^K, and Bq = 20/iG was calculated. 
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2.4 The CR proton and electron populations 

Fig. 3 sliows the volume integrated CR proton energy spectrum, E^-Np{E, t) = J 4TTr'^drE'^np{E, r, t), 



where E = mpC^ {y^T^Tf^ —1) , is the proton kinetic energy, and np{E^ r, t)dE = Anp^ f {p)dp = 
n{p, r, t)dp. We note that in the models with enhanced Bq {i.e., SI and S2) particles of charge 
Z can be accelerated to -Emax ~ Z ■ lO^^eV by th e early ST stage. This is consistent with 



Berezhko fc Volk 



20071) • In the early ST 



results of other analogous SNR simulations {e.g., 
stage the shock structure is significantly modified, with x = P2/P0 ~ 7 — 10, so E"^ ■ N{E) 
is concave upwards. This well-known behavior is a consequence of momentum dependent 
diffusion across the precursor; CRs near the injection momentum experience only the sub- 
shock velocity jump, while higher momentum CRs see a larger velocity jump. Thus, the CR 
spectrum is softer at lower momentum and harder just below the high energy cutoff than the 
strong shock test particle result, N{E) oc E'"^. Moreover, the Alfvenic drift in the precursor 
further softens the spectrum of newly accelerated CRs as the shock slows down and the 
Alfvenic Mach number decreases. 

It will be useful below to estimate the maximum CR proton momentum expected at a 
given time t/to- This we do by integr ating the standard express ion for the average momentum 



gained per unit time for a particle (ILagage fc CesarskyI 119831 ) 



dp 



Mo(l - Ma) -U2 



+ 



[u,{l-M-,')f 



p 



(7) 



1 + 



K2 90 



fto <7o— 3 ' 

where, once again, subscripts and 2 refer to upstream and downstream conditions in the 



shock frame, respectively. The useful parameter. 



(8) 



U0-VA-U2 x{^--Ma^)-i 

is the DSA test particle spectral index, including the upstream Alfven wave drift term, 
\uw\ = va = uq/Ma, for consistency with equation ([1]). The right-most expression for go 
neglects this term, since it is a relatively minor correction in our model SNRs. 

Neglecting the Alfven wave term and assuming a constant compression ratio, but 
including the diminishing shock speed in the ST solution, equation ([8]) gives for protons 
accelerated between tj and t, 



Pn 



0.54ulto 1 fto 



Kir. 



0.2 



0.2 



(9) 
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Figure 4. Volume integrated electron spectra in terms of particle kinetic energy at t/to = 1 (red dotted), 2 (green short 
dashed), 3 (blue long dashed), 6 (magenta dot-dashed), and 10 (black solid). (A colour version of this figure is available in the 
online journal.) 



For our simulations, wliich start from ti/to = 0.5, equation (|9]) asymptotes to Pmax ~ 
0.15(^0^2/^^) ~ 3,8 X 10^^{Eo/10^^y/^Bf,nJsS eV/c at large t/to using x = 4 for conve- 
nience. At t/to = 10 this gives, for example, Pmax ~ 7.8 x 10^^ eV/c for the SI model and 
Pmax ~ 5.8 X 10^^ eV/c for the S4 model. These estimates are consistent with the numerically 
determined cutoff energies in E"^ ■ Np{E) shown in Fig. 3. 

The CR electron spectra develop sub stantial differences from the proton spectra, despite 



their common PSA interactions (e.g. 



2009 



Blasi 



2010 



Webb, et al 



Zirakashvili fc Aharonian 



1984; 



Berezhko. Ksenofontov fc Volk 



20101 ). Fig. 4 shows the volume integrated CR 



electron energy spectra, E'^ ■ Ne{E,t) = An J r'^E'^ ■ ne{E,r,t)dr, at several times for each 
SNR model. Note that the electron spectra cutoff at much lower energies than the analogous 
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proton spectra shown in Fig. 3. The structures of the integrated electron spectra are complex 
around and below the cutoff. The low energy break in each spectrum identifies the energy at 
which the synchrotron/ Compton loss time for electrons downstream of the shock equals the 
SNR age. Below this energy radiative energy losses are negligible. The integrated spectrum 
between this break and the high energy cutoff is controlled by the downward movement of 
the electron energy cutoff as recently accelerated electrons move into the SNR interior away 
from the shock. The high energy cutoff itself develops at the shock and corresponds to the 
energy where synchrotron losses balance DSA gains. 

At the shock the electron cutoff momentum, pec, can be evaluated by balancing DSA 
and radiative losses. In particular, corresponds to the momentum at which the aver- 
age momentum gain by an electron in one pair of shock crossings, Ap/j^^ = (4/3)p[mo(1 — 
M.~]^) — U'2\/v, equals the momentum loss from radiation during the same period of time, 
Apr = h24:K2/[u2v\ + 6o4Ko/[(no(l - M']^)v\, where v = 
velocity, with 60 (p) and 62 (p) defined by equation ([3 
stream effective magnetic fields, respectively i^e.g., 
(test particle) result, 

(10) 



' TJif,/ {mpp)Y is the electron 
in terms of the u pstream and down- 



Webb. et all 



1984J ). This leads to the 



rrie 
rrir 



efffi , 



/ me 



1/2 



V o-tcB^ 



1/2 



R2 

X _|_ 1 go 

X ^'iff,o 90-3 



90 



.1 + ^-,0-3 \ 

where Us{t) = uq is the instantaneous shock speed. Here we assume kq/k2 = -B2/-B0 = X 
and the second expression corresponds to the limit, -Be// ~ B, which is applicable in our 



simulations. This Pec translates into a cutoff energy for CR electrons at the shock. 



3.3 PeV 



Bo,Mqo(l + X 



go ^ 
qo-3> 



Us(t) 



104 km s- 



(11) 



where the subscript /i on B indicates field strength in microGauss. 

Equation ( ITT]) evaluated with Bq ^j^ = 30 and x = 4 predicts E^c ~ 100 TeV for these 
models at time t ~ to, which is consistent with the simulation results shown in Fig. 4. This 
is also consistent with the upper limits on electron sp ectra derived from observation for a 



number of observed SNRs flRevnolds fc Keohane 



example, of 



I999I ) and with SNR simulation results, for 



Berezhko. Ksenofontov fc Volkl ( 12002| ). During the early ST stage, the maximum 



acceleration momentum gradually saturates as given in equation (|9]), but then it eventually 
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asymptotes to determined by the instantaneous shock speed as given in equation (ITT!) . 
As a resuh, for larger values of Bq, the electron cutoff energy is smaller, but the maximum 
proton energy is larger. 

Electrons advected into the SNR interior will continue to lose radiative energy with the 
combined synchrotron/IC radiative cooling time, tr{p) = p/b{p), given by 



Setting the cooling time behind the shock to the SNR age, tr(Pe,i) = t, provides a conserva- 
tive, rough estimate for the minimum electron momentum influenced by radiative cooling. 



Recall that Eg ^ p GeV. Below peb the integrated electron spectrum should be very similar 
in form to the proton spectrum at the same momentum, consistent with results shown in 
Fig. 3 and Fig. 4. 

To help interpret the integrated electron spectra, Ne{E), between the break momentum 
Pel, and the cutoff momentum Pec? we note in the thin shell approximation that Nf,{E) oc 
ne{E)r'^ Ar , and that the local electron spectrum, Ueip) = K^/pUplp), below p^t- In this case 
it is easy to show that the integrated electron energy spectrum below E^c would be steeper 
by about one in the spectral index than the integrated proton spectrum in the same energy 
range. The steepened slope would extend down to the energy Egb ^ peb GeV. On the other 
hand, in spherical flows following roughly the ST behavior, the postshock flow accelerates 
away from the high magnetic field region just inside the shock while reducing the magnetic 
field strength. These effects reduce the rate of cooling from that implied in equation [131 thus 
shortening the interval between p^c and Peb- Then the spectral steepening of the integrated 
electron population compared to the proton spectrum will be greater than unity, as observed 
in our results. 

3 MODELING OF CONTINUUM EMISSIONS FROM THE SNR MODELS 
3.1 The Cosmicp code for radiative processes 

The nonthermal radio to 7-ray emissions expected from CR electrons and from proton sec- 
ondary products in the simulated SNRs were computed through the post processing of model 
data using Cosmicp, a general purpose nonthermal continuum emissions code developed 
in-house around published radiative process and inelastic particle scattering formulations. 




(12) 




(13) 
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Cosmicp computes direct emissions from input electron/positron populations and calculates 
the inelastic collision products for photon and proton interactions with matter, including 
pion and lepton secondaries. Lepton secondaries are incorporated into Cosmicp, but they 
do not make significant contributions to SNR emissions, so are not included in our analysis 
here. Cosmicp was designed to be very flexible, so makes as few assumptions as possible. 
For example, the input energy distributions of energetic particles are arbitrary. Using inputs 
n{E,r) for CR electrons and protons along with the relevant environmental information, 
such as magnetic field, radiation field, and gas density, Cosmicp calculates the volume emis- 
sivity per unit frequency, j^{v,r), for each radiative process, where v is the emitted photon 
frequency. Bremsstrahlung from the thermal electron population is also included in our cal- 
culations with Cosmicp, but electronic emissions involving discrete atomic transitions are 
excluded. All calculations are in cgs units, which allows easy linking to actual observables. 



Specific CR electronic emissions incorporated into our analysis here include synchrotron, 
inverse Compton (IC), and bremsstrahlung processe s. The synchrotron , IC, and relativistic 
bremsstrahlung formalisms followed presentations in ISchlickeiserl (120021 ). The nonrelativistic 



bremsstr ahlung 
followed 



brmalism used to compute the spectrum due to the thermal population 
JacksonI (119991 ). The IC formalism uses an arbitrary, user supplied spectral form 
for the ambient photon field, so Cosmicp can model interactions with multiple blackbody or 
non-blackbody sources. It does, however, currently assume an isotropic incident radiation 
field. 



Hadronic interactions include, as noted, inelastic proton-proton and photon-proton inter- 
actions. Photopion produc tion in the latter case can be calculated from an arbitrary photon 
field using the formalism in lKelner fc AharonianI (120081). Pion generation from proton-proton 



interactions uses relations presented in 



Kelner. Aharonian. fc Bugayovl ( 120061 ) . The principal 



radiative process that results from hadronic interactions in our SNR context is secondary 
7r° decay, which is dominated by the 27 channel. Cosmicp can include helium in the hadron- 
hadron production of photons. However, that correction was not applied in the present calcu- 
lations. For a CR composition similar to that incident on the earth and an ISM with normal 
metalicity this correction would increase the pion production rate by ~ 25% Jschhckeiseil 



2OO2I). 
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3.2 Nonthermal emissions: Local interaction conditions 

Rates for nonthermal electron bremsstrahlung and p-p collisions leading to pion decay 7's 
depend on the product of the local plasma density and the CR electron and proton density, 
respectively. These are taken in our emission calculations directly from the simulations. 
Synchrotron and IC emissions depend on the spatial magnetic field distribution and the 
ambient radiation field respectively. As discussed in §2.2, we assume the local magnetic 
field strength scales with thermal plasma density; i.e., B{r,t) = BQ[p{r,t)/ pq]. The ambient 
photon field for our calculatio ns is uniforrn and models the Galactic radiation field in the 



solar neighborhood as given by ISchlickeiserl ( 120021 ). It consists of a sum of four components: 



the Cosmic Microwave Background (CMB) blackbody at 2.7 K with an energy density of 
0.25 eV/cm^, dust at 20 K with an energy density of 0.4 eV/cm^, old yellow stars at 5000 
K with an energy density of 0.3 eV/cm^, and young blue stars at 20, 000 K with an energy 
density of 0.09 eV/cm^. The total energy density for the ambient photon field is 1.04 eV/cm^. 
In evaluating equation ( IT2l) for electron energy losses this translates to an effective magnetic 
field of Br = 6.5 pG. We point out in §3.3.2, however, that IC scattering by the highest 
energy electrons of all but the CMB and dust emission radiation fields is limited by electron 
recoil, and reduced by the Klein-Nishina cutoff, so this value overestimates the effective B^. 
In our models electron losses are dominated by the synchrotron process, so this is a relatively 
small correction, in any case. 

3.3 Volume integrated radiation spectrum from the model remnants 

A general impression of the radiative emissions produced in the simulated SN la remnants 
can be obtained from Fig. 5, which shows for the SI simulation at t/to = 3 the volume 

requencies to PeV 7-rays. 
There is a narrow 7-ray 
window, 10^'^Hz < z/ < lO^^Hz dominated by bremsstrahlung from CR electrons interacting 
with thermal plasma. Above ~ lO^^Hz {E^ ^ 400keV) the spectrum is dominated by photons 
from 7r° decays for the ISM density assumed in the simulation. 

The spatially integrated radiative emissions at selected times for all four models are 
shown in Fig. 6. Independent of time, the radiation is mostly dominated by the electronic 
CR emissions at lower energies, but, for the most part, proton vr'' decay at high energies. 



integrated Spectral Energy Distribution (SED), uLi^ from radio 
Synchrotron emission dominates the spectrum for u < 10^°hJ^ 



2 Useful translation factors are: 1 keV 2.42 X lO^'^Hz, lO^'^Hz 0.41 keV. 
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Spectral Energy Distribution for SI at t/t^ = 3 

— I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — I — 




log[i. (Hz)] 

Figure 5. The volume integrated SED for the SI model at t/to = 3: total emission (solid black line), synchrotron emission 
(red dotted line), IC emission (green short dashed lino) bremsstrahlung (blue long dashed line), 7s from tt" decays (magenta 
dot-dashed line). (A colour version of this figure is available in the online journal.) 

The additional bump in the UV and X-ray seen at late times especially in model SI is due 
to thermal bremsstrahlung. This feature appears as the high energy end of the electron 
distribution is depleted by radiative losses, revealing thermal emission from hot, postshock 
plasma. The lower gas density in the S4 model reduces this thermal contribution, and also 
the nonthermal bremsstrahlung 7-rays. 

In the S4 model the 7-ray portion of the spectrum is a composite, dominated by 7r° 
decays between u ~ 10^^ - 10^^ Hz {E^ ~ 40 MeV - 400 MeV) and again above about 
V ~ 10^^ Hz [E^ ~ 40 TeV). However, the lower ambient gas density of the S4 model also 
reduces the vr'^ production enough that IC emissions are dominant between the synchrotron 
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Figure 6. The total volume integrated SEDs for models S1-S4 over time. The dotted red line is t/to = 1, the short dashed 
green line is 2, the long dashed blue line is 3, the dot-dashed magenta line is 6, and the solid black line is 10. (A colour version 
of this figure is available in the online journal.) 



cutoff and ~ lO^^Hz. IC emission is also predominant between i> ~ 10^^ — 10^'' Hz {E^ ~ 4 
GeV- 4 TeV). 

We outline below some useful relations for interpretation of the calculated radiative 
emissions in the context of these models, including some simple scaling relations that can 
be helpful in extending our results to different model parameter choices than we have used. 
We begin with discussions of the electronic emissions and follow with the essential elements 
to understand the calculated emissions produced by inelastic p-p collisions. 
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3.3.1 Electron Synchrotron Emission 

The electron synchrotron emission spans a frequency range from radio to roughly the emis- 
sion peak for electrons near the Lorentz factor 7ec ~ {mp/me)Pec- Above this energy the 
electron population is heavily modified by losses and begins to drop off severely. The asso- 
ciated synchrotron peak frequency is i^^c ~ 0.3 Vcriti where 



^crit 



'-UBlL^'ai{d) = 4.2 X 10^6 Hz 



(14) 



and ub is the nonrelativistic electron cyclotron frequency (IBlumenthal fc Gouldlll970l ). For 
the calculations presented here we set the angle of the magnetic field with respect to the 
observer, 6 = 90°. 

Using equations ffTTl) for the upper electron energy cutoff and f|T^ we can derive the 
following estimate for the cutoff frequency in the synchrotron spectrum in the limit that 



Beff ^ B] namely, 
^ 1.8 X 10^° Hz 



1 



go(l + X 



go ' 
90-3' 



10^ km s-i 



(15) 



For X = 4 and go = 4 the factor go[l + X'?o/(go -3)] = 17, so z/,c ~ 3 x 10^ ^ Hz fuJlO^km/s) 
consistent with the resu lts in Fig. 6. As noted by previous authors {e.g., 



Berezhko &: Vol 



5l3 



Bamba et al 



2003 



20041 ) the cutoff frequency is determined mainly the instantaneous shock 



speed and independent of magnetic field strength except through its implicit dependence on 
the factor x- This explains the similarity between the values for the synchrotron cutoffs for 
SI and S3 even though they have very different magnetic fields. 

The break frequency that corresponds to the estimated break momentum p^b depends 
only on magnetic field strength and the shock age as 

.3 -....\ 2 / n \ -3 



iy,b ~ 7.52 X lO'^Hz 



lO'^ yrs 
t 



B 



cfT 



100 /iG 



(16) 



The values of and Usb generated from equation ( |T5l) and (|T6l) . respectively, are con- 
sistent with the SEDs shown in Figs. 5 and 6. The photon energies corresponding to the 
synchrotron cutoff frequency fall in. the r ange 0.1 — 10 keV, consistent with observed Type 
la remnants f lReynolds fc Keohandll999l ). 

The synchrotron luminosities of all four simulations below the X-ray cutoffs can be 
approximately related to each other by a simple scaling relation that, in addition, pro- 
vides a means to extend these results approximately to different parameter choices. For 
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E <l{fi GeV the volume integrated electron energy spectrum shown in Fig. 4 can each be 
represented by 

N,{E,t)^{s-m,Et-^E-\ (17) 

where Ei ~ kBT2 represents the suprathermal injection energy of CR electrons, Ne is the 
volume integrated number of CR electrons and s is the energy distribution mean power law 
slope between Ei and < E < 10^ Gev. For reference, recall that momentum and energy 
distribution power law indices are r elated as, g = q — 2. The s ynchrotron luminosity can 



then be expressed approximately as (IBlumenthal fc Gouldlll970l ) 



Since Ne oc ^.K^/pMs^ and the injection faction ^ is nearly a constant for t/to > 1 in these 
simulations, Ng oc Msw {t/to)^^^. The injection energy, Ei scales roughly with the shock 
speed, so also as E^ oc Msw^"^ ■ In §2.2, we established the approximate scaling, x ^ Msw^^ , so 
B oc Msw^^ ■ Then uL^sync oc MsU''^'^^^^'^ in these simulations. For example, uL^sync oc mIw 
for s = 2 or vL^^gync oc mIw^ for s = 7/3. Since the index s varies slowly within the range 
9/5^s<7/3in these simulations, the radio luminosity should depend very weakly on Msw 
or t/toi consistent with the results in Fig. 6. 

Note that equation (fT8|) alone is not sufficient to estimate synchrotron emission behaviors 
within the X-ray band. The synchrotron emission in this band is generated predominantly 
by electrons near the spectral break Vgb and cutoff Vgc- Consequently, X-ray synchrotron 
emissions are heavily influenced by the positions of these frequencies. 



3.3.2 Electron Bremsstrahlung and IC Emission 

The bremsstrahlung and IC contributions are mostly subdominant in our models. At early 
times synchrotron emission dominates thermal bremsstrahlung over the bands where they 
both contribute. The bremsstrahlung peak is roughly Vnp ~ /CT2//1, while the synchrotron 
cutoff, Vsci is given by equation f|T5l) . Both cutoffs depend primarily on the shock speed 
squared, so evolve more or less together. The decreasing role of synchrotron emission in 
the X-ray band results, on the other hand, from the fact that CR electron energy losses 
downstream of the shock cause that population to become further depleted, dropping the 
relative emitted power between i/^c and Vsbi revealing the thermal bremsstrahlung emission. 



24 Edmon, Kang, Jones, and Ma 

Consequently for t/to > 6 in SI thermal bremsstrahlung is an important X-ray contributor 
(see Fig. 6). 

In another exception to the general emission behaviors, we note that IC emissions near 
/iz/ ~ 1 TeV in model S4 exceed the otherwise dominant pion decay gamma-rays. The S4 ttq 
decay emissions, on the other hand, are smaller by two orders of magnitude, reflecting the 
similarly lower S4 ambient density (see equation [24]). 

Some comments on the contributions and form of the IC SED may be useful. As noted 
in §3.2 the IC emissions result from a combination of incident photon fields with blackbody 
forms at different characteristic temperatures. The effective incident field is just their sum. 
While the CMB is a blackbody in both spectral form and energy density, the radiation from 
dust and stars has greatly diluted energy density compared to the Planck functions for their 
respective temperatures. In each case we can assign a frequency-independent dilution factor, 
-D, given by the ratio of the local energy density to the appropriate Planck function. The 
dilution factors for the radiation field properties listed in §3.2 are: for the CMB, Dcmb = 1, 
for thermal dust emission, = 5.3 x 10~^, yellow stars, Dy^, = 1.0 x 10~^^, and for blue 
stars, Dh^ = 1.2 x 10"^^ 

The IC scattering source function in the limit of Thomson scattering of incident photons 
with a Planck spectrum with color temperature, Tj, and dilution factor, Di, is 

TT ill C \ TT}jQ I 

w here rp is the classical electro n radius, -F(s) is a dimensionless scaling function as defined 



m 



Blumenthal fc Gouldl ( 1l970l ): for s ~ 2 -F(s) ~ 5. We have used the form of the electron 



spectrum given in equation ([Tj 

For a representative CR electron slope, s = 2, the different background radiation com- 
ponents contribute in equation 0191) according to factors oc T/^^Dj, which leads to a ranked 
list CMB, d, y*, b* represented in the ratios 1 : 0.58 : 0.03 : 0.004. This suggests that the 
CMB and the dust IR radiation are the predominant contributors to the IC SED, despite 
the roughly comparable photon energy densities of the stellar components (see §3.2). In this 
regard it is also worth noticing that the IC spectral peak at a few TeV involves photons 
scattered to energies close to E^^, the maximum possible, whereas the Thomson scattering 
limit Ej ~ ksTi^l^ would take photons into the PeV range when scattered from the incident 
stellar photon fields. In fact, photons from those higher temperature incident fields are scat- 
tered in the Klein-Nishina limit by the most energetic electrons, as given by the condition 
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on the incident photon energy, hu > rrieC^ /'-fee- This condition can be expressed practically 



as 



T, > i!^^ ^ 3.0 X 10^ (^) Kelvin. (20) 

For ~ 100 TeV only the CMB and, marginally, the dust radiation escape this limit. 
The scattering cross section in the Klein-Nishina limit is reduced from the Thomson cross 
section, so that, in fact, gamma ray IC emissions in these models are dominated by the CMB 
and dust emissions. 

We can use equation ( |T9l) to construct an approximate, simple scaling relation for the 
IC luminosity, vLyjCi in the Thomson regime assuming an electron distribution of the form 
in equation f|T71) . 

uL^.c oc ^eEr'utT'"^-'^'-'^'' (21) 

where the radiation energy density, Urad-, would be approximated by the sum of the CMB 
and dust radiation densities. The IC luminosity, uL^jc, scales as Msw similar to the 
synchrotron luminosity . In fact, the ratio of the two for emissions at synchrotron frequency 
^sync lec^B, and IC frequency, ujc <C lec^T^ in, say, the radio and gamma ray bands, is 



just 

i,T „{^+5)/2 /„ \(^-3)/2 

yL^,sync U%^'^" V ^IC ) ^ ^ 

This indicates that the ratio between the two contributions to the SED for each of our SNR 
models remains roughly constant over time. Since Urad is model independent, the ratio of 
IC gamma ray luminosity to radio synchrotron luminosity varies in these models roughly as 

oc l/Bl'\ 

3.3.3 7r° Decay Gamma-Ray Radiation 

Gamma ray decay products from inelastic p-p collisions, and ttq — 27, in particular, are 
commonly viewed as the "smoking gun" for hadronic DSA in SNRs. As noted above these 
emissions mostly dominate the gamma ray SEDs in the SNR models being discussed here. 
Although the detailed modeling of tt" decay gamma rays is fairly complicated, a simple 



analytic approximation offered by lPfrommer &: Enfilinl ( l2003l ) provides a good understanding 
of the 7r° gamma ray SED in these models. In particular, if the local CR proton energy 
spectrum above the threshold for pion production [E > 780MeV) is a power law, np{E) = 
(s - l)npE^-^E-', then the 7-ray source function, g^, (photons/time/volume/energy) is 
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n{s- l)^r % (23) 



where app ~ 3.2 x (0.96 + e4.4-2.4s-j ^ 10~^^ cm^ is an effective p-p cross-section for pion 
production, 5 = 0.14s~'^'^ + 0.44, while n is the density of the thermal protons, and = 
135MeV/c^ is the vr^ mass. 

This source function peaks for = m^c^/2 = 67.5 MeV. At low energies it asymptotes 
to a power law q^^ oc E^^, while it asymptotes at high energies to a power-law oc E~^. 
The source function cuts off above roughly 10% the proton high energy cutoff (Kelner et al. 
2006). The radiated power per unit volume (the SED) will scale with E^ g^. Using equation 
([23D for s > 2 the SED would peak around E^ = 67.5 MeV [{5 + 26/s)/{5 - 25/s)]i/(25)_ 
There is no SED maximum for s ^ 2, except at the high energy cutoff. 

Assuming a proton spectrum of the form in equation f|T71) one obtains a simple high 
energy scaling relation for the volume integrated pion decay luminosity 

i/L,,,o oc nNp E^-^ oc ^M,^E^~^ n E^~\ (24) 

This scaling is consistent with the results shown in Fig. 6, showing in particular why the 
pion decay emissions are similar in all the models except for S4 because of its lower rijsM- 

The proton distributions in these SNR models are not really true power laws, of course, 
but have concave upwards forms; that is, s decreases with energy. Protons near the pion 
production thresholds just above a GeV are mostly responsible for gamma rays near the low 
energy gamma ray peak. In that range typically s ~ 2.3 (see Figs. 3 and 4), The expected 
vr" gamma ray SED peak is near 1 GeV (z/ ~ 2 x 10^'^Hz), consistent with Figs. 5 and 6. 
The model vr^ SEDs in Figs. 5 and 6 also exhibit a second, higher energy peak. That results 
from the concave upwards character of the proton spectra and their eventual cutoffs around 
or above a PeV for t/to > 1. Except for model S4 the proton spectra have s < 2 below their 
cutoffs, making their upper vr*^ gamma ray peaks quite prominent. 

Since the the upper 7r° gamma ray peak is so prominent it can be used to determine the 



peak proton energy in the SNR {e.g., lAharonian et al.ll2007l ). From the gamma-ray produc- 



tion relations in Kelner et al. (2006) one can see that the peak in gamma-ray production for 
a monoenergetic proton of Ep > ITeV occurs at a photon energy E^ ^ O.lEp. Looking at 
Fig. 6 for S1-S3 the gamma-ray peak energy is around about 40 — 400 TeV (10^® — 10^^ Hz) 
which corresponds to a proton energy of 400 — 4000 TeV. This is consistent with the peaks 
in the proton spectra in Fig. 3. 
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Both IC and vr'' source functions depend on the number of CRs accelerated, while the 
pion-generated gamma rays also depend on the density of the thermal plasma. The ratio of 
Q'tt/q'/c found from equations ( l23l) and ( lT9l) . 



— TT^TT-, 7^ ^ Z - 



qic rl {kBTynieC^ 3s ym^ ) yksT 



s/S 

DF{s)K,,; ^ ' 

provides a useful comparison of the two processes for estimation and scaling purposes. 

Again, taking s = 2 as representative, then including CMB and dust contributions in 
the radiation field as discussed earlier, and setting K^/p = 10~^ along with n = 2 cm~^, 
appropriate to postshock conditions for model SI (see Table 1 and Fig. 2), we can use 
equation ( 125|) to estimate 

which is valid roughly for E^ ~ 400MeV {u ~ lO^'^Hz). This number is in reasonable 
agreement with the numerical results displayed in Fig. 5. Equation (j26l) suggests that 7r° 
decay emissions would need to be reduced by at least two orders of magnitude from those 
in model SI before IC emissions would become dominant in the GeV band. Note for scaling 
purposes from equation (!25l) that qn/qic oc n/ K^. ^. Consequently, GeV IC emissions would 
dominate when uism ^ 30i^e,p cm~'^, consistent with the model S4 results as displayed in 
Fig. 6. 



3.4 X-ray and gamma-ray brightness distributions 



Since X-ray telescopes and the new generations of gamma-ray telescopes have the ability to 
resolve images of SNRs with distances inside some tens of kpc, it is valuable to model spatial 
brightness profiles of these emissions. The profiles also provide insights to the similarities 
and differences among the different emission processes. The spatial profiles for all 4 models 
look very similar, so we will focus on model SI. The top panels of Fig. 7 show at several 



times the model S 



radial brightness p rofiles of broadband fluxes in the Chandra X-ray 



band (2 - 10 keV) ( IGarmire et al.ll2003l ) and the Hig h Energy Stereoscopic System (HESS) 



gamma-ray band (140 GeV - 40 TeV) (jBergd l2008l ). The profiles are plotted against the 
radius perpendicular to the line of sight from the center of the SNR, what we will call the 
projected radius i?, relative to the projected shock radius, Rs{t). The lower panels of Fig. 




Figure 7. The radial brightness profiles for the SI model in the Chandra X-ray band (2 - 10 keV) and HESS band (140 GeV 
- 40 TeV). Upper panels: Brightness profiles in terms of projected radius, R/Rs- Lower panels: The same profiles in angular 
units convolved to the resolution of Chandra and HESS with the SNR at 2 kpc. The line types arc the same as Fig. 5. (A colour 
version of this figure is available in the online journal.) 

7 present the same profiles in angular units when the SNR is placed at a distance of 2 kpc 
and the image has been convolved to the nominal resolutions in these bands for Chandra 
and HESS. 

Since the exponential cutoff of the synchrotron radiation spectrum, E-y{vsc) oc {t/to)~^^^ 
lies in the Chandra band, the 2-10 keV X-ray fluxes decrease quickly in time. This is evident 
in the evolution of the X-ray brightness profiles in Fig. 7. 

The X-ray brightness profiles are also very narrowly peaked at the projected shock po- 
sition. The X-ray profiles in Fig. 7 have characteristic apparent FWHM (Full Width Half 
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Max) of ~ 1%. This is consistent with X-ray observations of several young SNRs that have 
revealed thin nonthermal filaments at the forward shock with a FWHM thickness ^ 1% of 



the shock radius in projection (IBamba et al. 



2003 



Ballet 



2006 



Parizot et al 



20061). One 



can estimate the observed thickness of the shock rim when radiative losses are the domi- 
nant factor by the following relation from 



Parizot et al 



( 120061 ) ■ /^Rohs = 4.6u2tr, where t^ 



is defined in equation ( !T2|) . and the factor 4.6 is a geometrical correction that comes from 



projecting a spherical shock with an 



exponent ia 



the FWHM of the resulting profile ( iBalletl l2006l ). Taking U2 = Ust/x^ the ratio of shock 



emission profile onto a plane and taking 



speed to shock radius for the ST similarity solution, f/sT/''"ST = (2/5) /t, the relative FWHM 
shock thickness can be approximated as 

-2 / \ -1 / , \ -1 

(27) 



obs 



1.84 ftr 



X 



t 



0.11 



( gpff Y 

X \20012GJ 



le 



t 



SOOyrs 



2 X 10^^ 

where the numerical evaluation of tr from equation f|T2|) has been applied on the RHS. 
Taking x ~ 7 fo^' young SNR along with characteristic values for the rest of parameters, 
equation (j27|) indeed predicts ARobs/Rs ~ 1%- We note, in addition, applying equation [TT] 
in e quation ITTl that AR^h.^ oc indepe ndent of time or the SNR shock speed ((see, also. 



e.g. 



Volk. Berezhko. fc Ksenofontov 



2OO5I). 



The gamma-ray spatial profile in Fig. 7 is dominated by 7r° decay and is very broadly 
peaked at the shock. The peak at the shock is due to the emissivity scaling as the product of 
the CR density and the thermal gas density. Roughly, this will resemble j;^^^ oc p^, once again 
similar to thermal bremsstrahlung. However, the brightness profile is considerably broader 
than the emissivity owing to same projection effects noted earlier for the X-ray profile. 

The IC spatial profile, which is mostly buried under the neutral pion radiation in these 
models, will have a similar profile to the X-ray synchrotron, because they both come from the 
same population of electrons. Thus the same argument for why the synchrotron rim should 
be thin would apply to the IC spatial profile. However, the X-ray synchrotron profile will 
be slightly thinner than the IC profile owing to its dependence on the magnetic field. Since 
the IC rim will be thin and the vr*^ decay rim will be thick one could in theory differentiate 
between the two scenarios by looking at the shock thickness. If the shock is thinner in the 
TeV gamma-rays than the thermal gas distribution, like the X-ray synchrotron profile, then 
IC is dominant. But if the 7-ray profile has a thickness similar to the gas profile thickness, 
then 7r° decays are dominant. This level of discrimination, unfortunately, will not be available 
in the near future. 
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The lower panels of Fig. 7 provide a look at the SI model SNR placed at a distance 
of 2 kpc, about the distance to SN 1006. The intensity profiles have been convolved, as 
outlined below, to the angular resoluti on of the Chandra A CIS (Advanced CCD Imaging 



Spectromet er) in the X-r ay (0.5 arcsec 



Garmire et al. 



(120031)), and HESS around one TeV 



(3.6 arcmin iBergd (l2008l )). The Chandra ACIS has a diffuse flux sensitivity (3(T detection) 
of 1.6 X 10~^counts/cm^/sec/arcsec^ or about 1.3 x 10~^counts/cm^/sec/beam for 10 ksecs 
of observing time. Our model X-ray rims would be fully resolved by Chandra at the 2 kpc 
distance, although none of the earliest profile shown would be detected with a 10 ksec ACIS 
exposure. 

HESS resolution is sufficient to identify a shell-like structure of the gamma-ray profile, but 
would not be able to resolve the width of the rim, since only a few beams would fit across the 
source. The HESS diffuse flux sensitivity (8cr detection) of about 2.3x 10~^'^counts/cm^/sec/arcmin^ 



or ab out 9.3x10 counts/cm /sec/beam for 103 hours fl Garmire et al. 



2003 



Berezhko. Ksenofontov fc Vol 



20091 ) would be sufficient to detect this model source at each time shown. Detailed compar- 



ison to models of the gamma-ray rim will require the higher resolution of next generation of 
Imagin g Atmospheric Chere nkov Telescope's (lACT) such as the Cerenkov Telescope Array 
(CTA) Jwagner et al.ll2009l ). 



3.5 Comparisons to real SNRs 



We round out this discussion by referencing the emissions from these model SNRs to some 
real Type la SNRs. In doing so we emphasize that the models S1-S4 were not intended to 
reproduce any particular SNR, but rather to represent properties of CR modified Type la 
SNRs evolving in typical environments during the ST stage. The following exercise aims 
not to establish precise physical properties of the selected real SNRs, but to illustrate how 
broadband observations, along with simple scaling relations, can be used in conjunction 
with such models for a qualitative assessment of the physical state of an observed remnant. 
The potential for such applications is growing greatly with each new generation of high 
energy observatory. For instance several SNRs have been observe d by lACTs a t TeV photon 



energies, now including the Type la SN 1006 (lAcero et al.ll2010l ) and Tycho ( lAcciari et al 



201ll ). Similarly, Fermi co 



GeV gamma ray band {e.g., 



lecting an impres sive catalog of SNRs detected in the 0.2 — 300 



AbdoetaL 



20101). Fermi has not reported detection of any Type 



la SNRs so far. However, the core collapse SNR, Cas A, which has been seen at a distance of 
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3.4 kpc, provides a useful benchmark. Fermz measured a flux, vF,, ~ 5 x 10 ~^^ ergs/cm^/sec 
at 12.2cr signiflcance after 396 days of all sky observation (lAbdo et al.ll2010l ). This translates 
to an isotropic luminosity, vL^ = 7 x 10'^^ ergs/s. Over an equivalent time Fermi would not 
detect our remnants at the distance of Gas A. However Fermi is scheduled to observe in all 
sky mode for at least 5 years, so it could marginally detect (~ 3cr) models S1-S3 at t/to ~ 10 
at the distance of Gas A. 

As noted in §3.3.3, the TeV gamma-ray spectra of these model SNRs extend to photon 
energies well above the upper limits for lAGT's such as HESS of ~ 40 TeV. The proposed 
High Altitude Water Gherenkov experiment (HAWG) would be able to reach above this 

d for photons above 1 TeV of vFy ^ 1 x 



limit and would have a Sex de tection thresho 
10"^^ ergs/cm^/sec for 1 year (IGonzalezI |2008| ) . At a distance of 3 kpc that translates to 



a luminosity of uL^, = lO'^^ ergs/s. This would allow HAWG to detect all of our models 
excluding S4 at that distance after 1 year. 

For the other three primary bands, radio. X-ray and TeV 7-ray, actual observations of 
Type la SNRs already exist. We use the scaling relations in the previous sections to compare 
our models to two well observed Type la SNRs; SN 1006 and Tycho. We caution again that 
the simulations were idealized, so we should not expect close matches. Still, the comparisons 
may offer simple insights into the dynamical states of the SNRs and their GR populations. 



3.5.1 SN 1006 

SN 1006 lies at a di stance of about 2.2 kpc and at that d istance had a radius of 9.6 pc as 



of 1993 (t=987 yrs) (IBerezhko. Ksenofontov fc Volk 



2OO9I ). Ambient density estimates vary, 



but th e recent proper-motion-based dynamical estimate, ujsm ~ 0.085 cm 



Katsuda et al. 



(120091 ) is representative. In order to establish a dynamical state for the comparison we utilize 
the ST expressions in §2.2 to obtain 

The age, t, is well constrained, and the shock radius, r^, is relatively well constrained 
by observation. Thus, even though po and are not individually as well constrained, their 
ratio po/ Eo is flxed at least within the ST paradigm. Since the explosion energy is likely close 
to lO^-'^erg (one "FOE"), by assuming an ejected mass, M^j = 1.4M0, we c an express the 



ambient density constraint as ujsm ~ 0.037 EpoE cm ^, consistent with the 



Katsuda et al 



( I2OO9I ) estimate. Similarly, to = 520 EpQ^ yrs, so at the observed time, t/to = 1-9 E' 



5/6 

FOE- 
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The swept up mass, Mg^ ^ A.QMqEpqe- We can use our simulation SI for comparison with 
SN 1006. 

The model SI 1.4 GHz radio luminosity, vLy ^ 3 x I O^^-Efqk ergs/sec. The ni easured 
1.4 GHz flux of SN1006 is vFy ^ 2.2 x lO'^^ ergs/cmVsec ( Reynolds fc EllisoJ 1992 ). which 
translates to an isotropic luminosity, vLy ^ 1.3 x 10"^^ ergs/sec, so slightly larger than the 
model. As noted in §3.3.1 the model radio emission is almost independent of time, or swept 
up mass, but it does scale with the electron injection efficiency, Kf,/p and the magnetic field, 
(assuming s = 2). So, a modest increase in this combination would bring the 
two luminosities into satisfactory agreement. 

The expected X-ray synchrotron luminosity is set within the model by the radio luminos- 
ity and by the break frequency, Ugb oc Eo/mIw . For model SI at tjto^l the br eak frequency. 



Bamba et al 



10 Hz, compared to the "roll-off" frequency in SN 1006 estimated by 
(120081 ) to be z/j-oz; ~ 6 x lO^^Hz. Fig. 6 shows the model SI X-ray luminosity near the break 
to be roughly lO^^erg/sec. The measured ( 3.1 — 2 keV (mostly nontherrn al) X-ray flux of SN 



1006 is vFy ^ 1.42 x 10"^" ergs/cmVsec (lAUen. Petre fc Gotthel 



200ll ). This translates to 



an isotropic luminosity, vLy ^ 8.1 x 10'^^ ergs/sec. The X-ray comparison is consistent with 
the radio comparison; the best match would come from a modest increase in K^j^B'^l'^ . 



HESS has reported detection of VHE gamma-ra ys from SN 



006 with a flux above 1 



TeV gamma ray flux, vF^, 3 x 10"^"^ ergs/cm^/sec (lAcero et al.ll2010l ). This translates to 
a luminosity of vLy ^ 2 x 10'^^ ergs/sec. This is very close to the (pion-decay) luminosity 
for model SI shown in Fig. 6. 



3.5.2 Tycho 



Tych o's SNR may not yet have reached the ST dynamical stage (jVolk. Berezhko. fc Ksenofontov 



20081 ). but we include it i n orde r to ill ustrate another application of the scaling relations. In 



this context we note that 



Kand (I200I found evolution of the outer, blast wave, shock to be 



represented reasonably well by ST scaling well before the actual ST phase begins. 

Tycho has a distance that is variously estimated, but mostly less than about 3.6 kpc. 



We will parameterize our comparison including a distance 



respo nding radius in 1983 (t=411 yrs) was = 4 d^^Q pc (IVolk. Berezhko. fc Ksenofontov 



'actor d^.H = d/3.6k'pc. The cor- 



20081 ). The SNR lies in the plane of the galaxy in a warm ISM region of estimated density 
~ 0.3 cm~^. The warm ISM context makes models SI - S3 the most appropriate compar- 
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ison. Since we have established the vahdity of the emission scahng relations among these 
models, we focus on model SI. We follow the same procedure as for SN 1006, except we 
now include distance as an adjustable parameter in order to illustrate its influence on the 
comparison. Equation (l28l) leads in this case to ujsm = 0.43 EpoE/d^.Q cm~^ along with 



t/to = 1.8 E]l^j,/dy^ and M,^ ^ A.^qEfoe/ dl 



Tycho's measured 1.4 GHz radio flux is z/F^ ^ 1.6x10 ergs/cm /sec ( iReynolds fc Ellison 



19921 ). which translates to an isotropic luminosity, uLy ^ 2.5 x 10'^^ (i|g ergs/sec. The scaling 



given in equation f|T8l) using nominal parameters leads to a model SI 1.4 GHz luminosity 
estimate appropriate for Tycho, uL^ ^ 4.9 x 10^^ EpoE/dl ^ ergs/sec, or about five times 
smaller than the observed luminosity. If we wanted to achieve a better match using the full 
scahng from equation (JT8l) . we would increase the factor ^Ke^pB^/'^EpoE/ d^^ by this factor 
of five. 

Additional information comes from the X-ray and 7-ray porti ons of the SEP. Tycho's 



meas ured 10 keV X-ray flux is pF^ ^ 3.2 x 10 ergs/cm /sec (lAllen. Gotthelf fc Petre 



I999I ). This translates to an isotropic luminosity, vLy ^ 4.9 x 10'^^(i|g ergs/sec. The model 
SI flux at 10 keV would be for nominal parameters, uL^/gync ~ 10^^ EpoE/d^^ erg/sec. 
This is once again about a factor five smaller, so leads to a consistent way to compare the 
model to the observations. The Very Energetic Radiation Ir naging Telescope Ar ray System 



(VERITAS) group recently reported a 7-ray flux above 1 TeV lAcciari et al.l (120111 ). uFy ^ 3x 
]^Q-i3 ergs/cm^ /sec. This translates to an isotropic luminosity, uLy f» 4.7 x 10^^ (i|gerg/sec, 
consistent with the properties of our model SI. 



4 SUMMARY 

Using kinetic, nonlinear DSA simulations, we have calculated the energy spectrum of CR 
protons and electrons accelerated at blast waves from Type la supernovae and the emission 
contributions of those cosmic rays during Sedov- Taylor phase. The DSA model is based on 
the assumptions that streaming instability amplified Alfven waves lead to Bohm-like diffu- 
sion, and that particle escape at the highest energy is negligible in the time frame of interest. 
With the CR proton injection rate from the thermal plasma, E, ~ 10"'^, obtained here, the 
strong SNR shocks become significantly modified by the formation of a CR precursor that 
pre-compresses the inflowing plasma. This both weakens the dissipative subshock and, com- 
bined with Bohm diffusion, leads to a CR proton spectrum that is steeper at low energies 
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than the strong shock test particle form, Np{E)dE oc E~'^dE, but flatter at high energies. 
In addition, the drift of scattering centers with respect to the bulk plasma also softens the 
CR spectrum and the dissipation of scattering waves reduces the CR acceleration efficiency. 
These effects become more important for stronger upstream magnetic fields and for a lower 
ISM density, because the Alfven speed is va oc B/^/nisM- 

We assume the electron to proton ratio is K^jp = 10~^ that is, ne(pinj) = -f^e/p'^plPinj) 
near the injection momentum. Although the value of Ke/p ~ 10^^ is implied by the Galactic 



CR flux observed at Earth, the smaller value is preferred 



mult 


i- 


2009 





jand photon spectra of severa . 



Morlino. Amato fc Blasi 



2009 



young SNRs (e.g., 



'or source spectra that best flt 



Berezhko. Ksenofontov &: Volk 



Zirakashvili fc Aharonian 



20101 ). However, since the 



electrons are passive with respect to the hydrodynamics the results of this paper can be 
rescaled to any (small) K^/p. 

The thermal and nonthermal continuum emissions from the model SNRs were computed 
using a radiative process code Cosmicp. The spectra are dominated in most cases by CR 
electron emissions in the X-ray band and below. Gamma-ray emissions were dominated in 
most cases hy p — p produced 7r° decay. For the hot, low density ISM case inverse Compton 
electronic emissions dominated the TeV 7-ray band, however. We note of course that the flnal 
comparison of hadronic versus leptonic origin of GeV TeV 7-rays depends on the ambient 
matter density and the electron to proton ratio as Uism/ Ke/p aniong other parameters. 

We presented simple scaling relations for the synchrotron, IC and vr^ decay luminosities 
that seem to provide appropriate approximate means to rescale the numerical results to 
different times and model parameters. We applied these relations as simple tests to SN 1006 
and Tycho's SNR demonstrating a basic consistency and a means to explore parameter space 
dependencies. 
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